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Abstract. In this letter, we p resent numerical calculations made to investigate the possible resolution dependence 
of the Shakura & Sunyaev ( 1973 ) viscosity parameter a from local magnetohydrodynamic simulations of the 
magnetorotational instability (MRI) . We find that the values of a do indeed depend significantly on the numerical 
resolution but also that when the highest resolutions attainable by the computational resources available are used, 
the growth of the a-parameter seems to saturate. The values of a are at most of the order of 10 -3 , which indicates 
that the sole presence of turbulence due to dynamo generated magnetic field in the di sk is n ot enough to reproduce 
as of the order unity which could explain some observational results (e.g. Cannizzo 1993). 
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1. Introduction 

The biggest problem in the theory of accretion disks has 
been the lack of knowledge of a widely applicable physi- 
cal process which could reproduce the observed effective 
angular momentum transport. A decade ago, the magne- 
torotational instability was realised to have signifigance 



in the context of accretion disks (Balbus & Hawley 1991). 
The MRI was seen to excite and sustain turbulence giving 
rise to turbulent viscosity and outward angular momen- 
tu m transport in local MHD calc ulations (e.g . Hawley et 
al. |1995[ |1996| ; Brandenburg et al. |l995| , |l996| ; Stone et al. 
1996|). 



The efficiency of this process can be measured in 
terms of the Shakura & Sunyaev a-parameter which 
parametrises the turbulent viscosity as v t = ac s H, where 
c s is the sound speed, and H the height of the disk. The 
above parametrisation was made assuming the motion of 
the gas to be subsonic and that the turbulent eddies are 
smaller than the disk height, resulting in a as a dimen- 
sionlcss number whose values vary between zero and one. 
The interpretation of the lightcurves of some dwarf novae 
indicates that the value of a should be of the order 10~ 2 — 1 



(e.g. Cannizzo 1993) 



However, in the aforementioned numerical studies the 
values of a have been at least an order of magnitude too 
low to account for the observational results. On the other 
hand, the numerical resolution was seen to affect the val- 
uesof a (e.g. Brandenburg et al. 1996 ; Miller & Stone 
[2000 ). Doubling the resolution increased the value of a 
by a factor of ~ 1.5 in these studies. This possible de- 
pendence has not been studied systematically yet, casting 
some doubt that high enough as can be obtained from 
local MRI calculations. Therefore, the purpose of this pa- 
per is to perform the previously missing survey on the 
resolution dependence of the a-parameter. 

2. The MHD model 

The computational domain in our calculations is a rect- 
angular box with dimensions Li, i — x,y,z, at a radius 
i?o from the centre of rotation. The distance to the centre 
of force is large in comparison to the box dimensions (i.e. 
Rq Li), under which assumption new locally cartesian 
coordinates can be taken into use and the equations of mo- 
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tion can be linearised (e.g. Spitzer & Schwarzschild 1953). 
This simplification is usually referred to as the shearing 
box approximation, in which we solve a set of non-ideal 
MHD-equations in cartesian coordinates (see Caunt & 
Korpi |2001| ). 

We apply small molecular viscosities everywhere and 
additionally artificial shock- and hyperviscosities in order 
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Fig. 1. Turbulent kinetic and magnetic energies from the run Rl. The different components are denoted by solid 
(radial), dotted (azimuthal), and dashed (vertical) lines. The total energy is denoted by dash dotted line. 



Table 1. The values of maximum and time averaged a— 
parameters for the two sets of simulations. Time averages 
were taken up to 13 T or b (duration of the largest calcula- 
tion), where T or b = 27r/Slo is the orbital period. 



Initially the disk is in hydrostatic equilibrium which yields 
a Gaussian density stratification 



hip = lnpo - jp ■ 



(2) 



Run 


Resolution 


a max [10" 3 ] 


<a> t [KT 4 ] 


Tl 


15 x 15 x 31 


2 


3 


T2 


31 x 31 x 63 


13 


10 


T3 


47 x 47 x 95 


22 


16 


T4 


63 x 63 x 127 


23 


32 


T5 


95 x 95 x 191 


32 


39 


T6 


127 x 127 x 255 


28 


41 


Rl 


31 x 31 x 63 


14 


7 


R2 


47 x 47 x 95 


16 


18 


R3 


63 x 63 x 127 


32 


26 



The initial magnetic field configuration is a sinusoidally 
distributed vertical field in the radial direction giving a 
zero net field. In terms of the vector potential this config- 
uration can be represented as 



L x [2-kx\ „ 

A = — Bn cos v . 

2tt ° V L x J y 



(3) 



to resolve shocks and (unphysical) small scale oscillations, 
respectively. The artificial viscosities also serve the pur- 
pose of stabilising the numerics. 

We adopt periodic boundary conditions in the radial 
and azimuthal directions, and closed insulating stress free 
boundaries in the vertical direction. The radial boundary 
takes into account the basic shear flow. Numerical imple- 
mentation of the model described above is presented in 
the paper by Caunt & Korpi ( 2001 ). 

3. Physical setup 



dimensionless quantities as 
i.e. length is measured 



We adopt the same 
Brandenburg et al. ( 1995 ) 
in units of the initial density scale height [x] = Hq, time 
in units of [t] = (GM/Hq)^ 1 / 2 , and the density in units 
of [p] = po, the initial density at the midplane of the disk. 
The unit of magnetic field will then be [B] = [u^poPo) 1 ^ 2 , 
where [u] = [x]/[t\. Dimensionless quantities are obtained 
by prescribing 



H = GM = p a =p = l. 



(1) 



where Bq is calculated from the plasma beta, (3 = 
2pqp/Bq, which has an initial value of 100 in the mid- 
plane of the disk. Random velocity perturbations of the 
order 10 _3 c s were added in the velocity field to start up 
the instability. 

4. Calculations 



Two sets of calculations were made, differing in the box 
size. The numerical resolution was varied systematically, 
see Table |l| for complete details. In the R-set the used 
box size was the standard 1 x 2tt x 4 in units of the ini- 
t ial de nsity scale height Ho, as used b y e.g. Hawley et al. 
( 1995 ) and Brandenburg et al. ( 1995 ) to be able to com- 
pare the results with previous studies and thereby validate 
the code. In the T-set the smallest box size in which the 
MRI was still seen to be active and where the effects of the 
boundaries were not yet significant, namely 1x1x2, was 
chosen in order to achieve the best possible spatial resolu- 
tion. The lowest resolution runs were advanced up to ~ 150 
orbits, whereas the largest calculation could only be con- 
tinued until 13 orbits in a reasonable amount of time. The 
duration of the largest calculation determined the length 
of the time average for the energies and the a-parameter. 
In the calculations we systematically increase the resolu- 
tion keeping the box size and other parameters fixed, and 
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Fig. 2. Maximum and time averaged energies as a function of mesh size, Aj = Li/ni, where ni is the number of 
gridpoints and i = x,y, or z. The errorbars are given by a constant factor divided by the number of datapoints used 
in the calculation. 



calculate the resulting Shakura-Sunyaev a-parameter in 
the following two ways. Firstly, we calculate a using the 
M axwell and Reynolds stresses (see e.g. Brandenburg et 
al. |1995|) 



Ostrc 



< pU X Uy > 



-/!q < B X By > 



Secondly, we use the rate of dissipation 



a. <p e >)c 



(| n ) 2 c s0 Ho <p > 



(4) 



(5) 



In Eqs. (|) and (g) the square brackets denote volume 
averages, and we use the initial values of density, sound 
speed and density scale height as the normalisation factor. 
As the values of a from both methods coincide rather well, 
in the following we give only results from the first method. 
It is worth noting that the values of a calculated in this 
way are smaller than the values of e.g. Hawley et al. (1995 



1996) by a factor of 3/\/2 due to different normalisation 



factors. 



5. Results 



The overall behaviour of the system is clearly visible from 
Fig. ^ where we plot the kinetic and magnetic energies for 
the run Rl. The behaviour of the energies in the higher 



resolution simulations is very similar to the run Rl, ex- 
cept that the energies are growing with resolution, see 
Fig. H After the initial exponential growth, the energies 
saturate into a self-sustained turbulent state. Although 
large fluctuations are seen, no persistent decay or growth 
is seen after about five orbits. The initially weak mag- 
netic field is amplified considerably due to dynamo action 
which leads into a configuration where the field is strongly 
dominated by the azimuthal component, followed by sig- 
nificantly weaker radial and vertical components. The be- 
haviour of the magnetic energy is very similar to the ear- 
lier studies (e.g. Brandenburg et al. 1995). The kinetic 
energy, however, is slightly dominated by the radial com- 
ponent in our simulations, whereas in some earlier works 
the azimuthal component has also dominated the kinetic 
energy. Even though the magnetic energy is at all times 
larger than the kinetic energy, the system is dominated by 
thermal energy which is usually of the order of ~ 100 times 
larger than the magnetic energy. This causes the ratio of 
thermal to magnetic pressure to be such that the MRI is 
active at all times in the simulations. 

Other characteristic numbers describing the system 
include the ratio of the Maxwell and Reynolds stresses, 
which are the main transport terms extracting energy 
from the Keplerian flow and transforming it into turbulent 
energy. In previous studies this ratio has been between 3 
and 6. The value found in our simulations is close to 3, 
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Fig. 3. Time averages of the a-parameter as a function of mesh size, Aj, where i = x,y, or z. 



being in satisfactory agreement with the previous results. 
Although the results are not exactly identical to the earlier 
ones, the differences are not huge, and the basic character- 
istics can be well described by our model, indicating that 
it is valid. 

The maximum and time averaged energies as functions 
of the mesh size, Aj = Li/ni, for the T-set are shown in 
Fig.^|. The maximum and time averaged kinetic energies 
both seem to continue almost linear growth beyond the 
resolutions investigated here, whereas the magnetic ener- 
gies seem to saturate at the highest resolutions. A very 
similar trend was also found in the R-set (not shown), 
which indicates that the small box calculations are a valid 
tool in the study of resolution dependencies. 

In Fig.|3] the maximum and time averaged a— 
parameters as functions of mesh size are shown for the 
T-set of runs. The actual numerical values of a are given 
m Table § As can be seen from this table for the R-set, 
doubling the resolution (from Rl to R3), the value of a 
increases by a factor of ~3, which is even more dram atic 
than in previous studies (e.g. Brandenburg et al. 1996|) . In 
the T-set, the strong growth seen in the lowest resolutions 
seems to stop when resolution is increased further from 
the run T4. For the largest calculations, a is almost con- 
stant. The largest values of time averaged a we obtained 
were of the order of 0.004 which are roughly two orders 
of magnitude too low to account for the aforementioned 
observational results of dwarf novae. 

6. Conclusions 

In this letter, we have numerically investigated the resolu- 
tion dependence of the Shakura & Sunyaev a-parameter 
from a local MHD model of a weakly magnetised accretion 
disk where the MRI is acting. We find that the numerical 
resolution affects the energies and the values of a signifi- 
cantly at lower resolutions. However, the strong growth of 
magnetic energy stops at the highest resolutions indicat- 
ing that the Maxwell stress cannot extract more energy 
from the Keplerian shear flow when the resolution is in- 
creased which is also reflected in saturating trend seen in 
the values of a. As the largest values of a obtained in this 
study were only of the order 10~ 3 , the obvious conclusion 



is that it is not possible to obtain high enough as due to 
the MRI, at least not from local restricted models as the 
one studied here. There are, however, recent global stud- 
ies of the MRI (e.g. Hawley 2001 ), which models exhibit 
large scale wave phenomena producing significantly higher 
values of a (of the order ~0.1). 

Another hydromagnetic instability, namely the accre- 
tion ejection instability (AEI), investigated by Caunt &; 
Tagger (2001) has recently been identified as a source of 
angular momentum transport, producing similar physical 
characteristics (spiral waves and comparable values of a) 
as the global MRI simulations, but which instability can- 
not operate without a fairly strong magnetic field. The 
AEI cannot generate this field by itself so it has to orig- 
inate either from an external source or be generated by 
some intrinsic dynamo process, such as the one caused 
by the MRI. This suggests that the MRI is still a viable 
candidate to be responsible for the angular momentum 
transport (possibly indirectly via the AEI) in accretion 
disks and that it is the restricted geometry of the local 
models behind the too small values of a. 
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